%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Explanation of what the file does
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% This function obtains a new guess for the path of utility dots in the RUV
% model with forward looking agents
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Inputs needed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 1. Matrix of sectoral price dots (PHma), size I*(S-1)*T
% 2. Final consumption shares (Amat), size I*(S-1)
% 3. Matrix of sectoral wages dot (Wmat), size I*(S-1)*T
% 4. Matrix of sectoral actual labors dot (Lmat), size I*(S-1)*T
% 5. Matrix of sectoral labor supplies (EllEra), size I*S*T
% 6. Mobility matrices (MuEra), size (I*S)*(I*S)*T
% 7. Discount factor (ba), scalar
% 8. Migration elasticity (ka), scalar
% 9. Sectoral elasticity (nu), scalar
% 10. Initial guess for U dots (UInit), size I*S*T
% 11. Degree of risk sharing (zrs), scalar
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Outputs produced
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 1. New path of utility dots (UNew), size I*S*T

function UNew = NEBA_Block3(PHma,Amat,Wmat,Lmat,EllEra,MuEra,ba,ka,nu,UIni,zrs)

% Define some important variables

I                = size(EllEra,1);
S                = size(EllEra,2);
T                = size(EllEra,3);
EllRel           = EllEra(:,2:end,:);
EllHat           = EllRel(:,:,2:end)./EllRel(:,:,1:end-1);
Pagg             = repmat(prod(PHma .^ repmat(Amat,1,1,T-1),2),1,S-1,1);

% Things for imperfect risk sharing

LEra             = NaN(I,S-1,T);
LEra(:,:,1)      = EllRel(:,:,1);
LEra(:,:,2:end)  = cumprod(Lmat,3).*repmat(EllRel(:,:,1),1,1,T-1);
unempprob        = 1-LEra./EllRel;
Deltars          = zrs.^(-unempprob) .* ((1-unempprob*zrs)./(1-unempprob)).^(unempprob-1);
Delrsdot         = Deltars(:,:,2:end)./Deltars(:,:,1:end-1);

% Omegas and utilities

Omega            = ones(I,S,T-1);
Omega(:,2:end,:) = (Wmat.*Lmat)./(Pagg.*EllHat.*Delrsdot);
UNew             = ones(I,S,T);

% Run the backward induction loop where the future u dots are used to
% obtain the current u dots (the sectorial compensation omegas and mobility 
% matrices are also used)

for counter = 1:T-1
    t            = T-counter;
    Ufut         = UIni(:,:,t+1)';
    MuHash       = MuEra(:,:,t) * kron(eye(I),ones(S,1));
    MuCond       = MuEra(:,:,t) ./ kron(MuHash,ones(1,S));
    MuCond(isnan(MuCond)) = 0;
    Step1        = (repmat(Ufut(:)',I*S,1)).^(ba/nu) .* MuCond;
    Step2        = (Step1 * kron(eye(I),ones(S,1))).^(nu/ka);
    Step3        = sum(MuHash .* Step2,2).^ka;
    Step4        = reshape(Step3,S,I)';
    UNew(:,:,t)  = Omega(:,:,t) .* Step4;
end

end
